AD-A202  693 


erne  FILE  COW 


r£ 


☆ 


☆ 


W( 


DTIC 

^electe 

%,  JAN  1  81988 


I 


w  right-Patte  rson  Air  Force  Base,  Ohio 


DISTRIBUTION  STATEMENT  A 
Approved  far  public  releaee; 

T  T  » 1*  »r>  • 


o  O 


AFIT/GA/AA/88D-7 


OPTIMAL  LAUNCH  TIME  FOR  A  DISCONTINUOUS 
LOW  THRUST  ORBIT  TRANSFER 
THESIS 

Jeffrey  M.  McCann 
Captain,  USAF 

AFIT/GA/AA/88D-7 

DTIC 

s 

H 

Approved  for  public  release;  distribution  unlimited 


AFIT/GA/AA/88D-7 


OPTIMAL  LAUNCH  TIME  FOR  A  DISCONTINUOUS 
LOW  THRUST  ORBIT  TRANSFER 

THESIS 

Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 

In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science  in  Astronautlcal  Engineering 


Jeffrey  M.  McCann,  B.S. 

Captain,  USAF 

December  1988 

Approved  for  public  release;  distribution  unlimited 


Eieiacg 


Electric  propulsion  is  a  viable  alternative  to  chemical 
propulsion  for  orbit  transfers.  The  difference  between 
present  day  chemical  rockets  and  the  near  term  electric 
rockets  is  similar  to  the  difference  between  a  Lamborghini 
Countach  and  a  moped.  Both  will  get  you  to  where  you  want  to 
go  but  the  moped  is  much  more  efficient  and  considerably 
slower.  Electric  rockets  could  be  the  moped  of  transfer 
vehicles . 

Electric  propulsion  will  improve  the  cost  efficiency  for 
payloads,  thereby  making  this  endeavor  more  profitable.  As 
budget  constraints  mount,  the  space  industry  must  look  for 
more  cost  efficient,  yet  safe,  systems  for  deploying  payloads 
to  various  orbits.  Electric  rockets  may  be  the  solution  to 
this  system  search.  However,  this  solution  hinges  on  one 
notion:  transfer  time  is  not  important.  The  slogan  for  a 
NASA/USAF/industry  electric  rocket  could  be:  if  you've  got 
the  time,  we've  got  the  rocket. 

I  would  like  to  thank  my  advisor  Dr.  William  wiesel  for 
Introducing  me  to  this  up  and  coming  technology  and  for  his 
patience  during  my  learning  curve.  I  would  also  like  to 
thank  my  committee  members  Capt  James  Planeaux  and  Capt 
Rodney  Bain.  A  special  thanks  is  due  to  Capt  Bain  for  his 
advise  and  guidance  while  I  was  a  part-time  AFIT  student 
prior  to  my  full-time  acceptance. 
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Abstract 

This  study  searches  for  the  optimal  launch  time  of  a 
discontinuous,  low  thrust  transfer  between  non-coplanar , 
circular  orbits.  The  spacecraft  is  assumed  to  be  a 
solar-powered  rocket  that  cannot  provide  thrust  in  the 
earth's  shadow. 

Two  timescales  are  used  to  calculate  a  minimum  time 
trajectory.  The  fast  timescale  produces  a  control  law  which 
maximizes  a  change  in  orbital  elements  for  a  single  orbit. 
The  slow  timescale  incorporates  the  control  law  so  that  the 
final  boundary  conditions  are  met  in  minimum  time.  Minimum 
transfer  times  are  determined  between  winter  and  summer. 


ix 


OPTIMAL  LAUNCH  TIME  FOR  A  DISCONTINUOUS, 

LOW  THRUST,  ORBIT  TRANSFER 

INTRODUCTION 

The  orbit  transfer  condsldered  in  this  study  Is  between  a 
200  km  circular  orbit  (a  %  1.03  DU)  with  an  inclination  of 
28.5  degrees  and  a  geosynchronous  orbit  (a  %  6.6  DU).  The 
launch  time  is  broken  into  two  components:  time  of  day 
(i.e.,  ascending  node)  and  time  of  year.  The  optimal  launch 
time  is  determined  by  minimizing  the  transfer  time. 
Discontinuous  low  thrust  is  possible  using  a  certain  type  of 
electric  propulsion.  An  overview  of  electric  propulsion  is 
given  below. 

Propulsion  systems  may  be  classified  into  two  basic 
categories : 

1.  Endogenous,  which  use  energy  stored  within  the 
propellants  to  create  thrust.  Solid  rockets  and  liquid 
rockets  are  common  examples  of  endogenous  systems. 

2.  Exogenous .  in  which  the  energy  is  supplied  to  the 
propellant  from  an  outside  power  source.  All  electric 
propulsion  systems  are  exogenous,  although  some  can  be  a 
combination  of  the  two. 

The  most  significant  advantage  of  an  exogenous  system  is 
that,  if  external  energy  is  available  for  accelerating  a 
propellant,  the  resulting  specific  impulse  can  be  much  more 
than  that  of  an  endogenous  system.  For  example,  an  ion 
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thruster  system  with  a  specific  Impulse  of  3000  seconds  would 
require  only  2000  kg  of  ion  propellant  as  compared  to  15,000 
kg  of  chemical  propellant  for  a  liquid  oxygen-liquid  hydrogen 
upper  stage  rocket  with  the  same  total  impulse.  The  dry 
weights  of  the  two  systems  are  also  similar,  resulting  in  a 
significant  advantage  for  the  ion  thruster  system  (7:xii). 

Electric  propulsion  devices  are  low  thrust  devices.  This 
allows  very  fragile  large  space  structures  to  be  transported 
from  low  earth  orbit  to  geosynchronous  orbit.  In  addition, 
since  the  propellant  is  a  small  fraction  of  the  overall 
system  mass,  weight  growth  of  the  payload  during  the 
construction  phase  of  the  project  can  easily  be  accommodated 
by  thrusting  for  a  longer  period  of  time  (7:xiil).  Payload 
mass  is  not  a  constraint  as  in  chemical  systems,  it  merely 
requires  longer  transfer  times.  Other  benefits  include 
better  cost  efficiency  (lower  dollar/kg  of  payload)  and  less 
complex  engine  design. 

There  are  three  classes  of  electric  propulsion  devices, 
all  of  which  are  capable  of  high  specific  Impulse. 

1.  Electromagnetic .  The  electromagnetic  thruster  is 
often  called  the  plasma  thruster.  In  this  thruster,  the 
propellant  gas  is  ionized  to  form  a  plasma,  which  is  then 
accelerated  rearward  by  electric  and  magnetic  fields. 

2.  Electrostatic .  The  best  known  type  of  electrostatic 
thruster  is  the  ion  thruster.  As  in  the  plasma  thruster, 
propellant  atoms  are  ionized  by  removing  an  electron  from 
each  atom.  In  the  electrostatic  thruster,  however,  the 
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electrons  are  removed  from  the  Ionization  region  at  the  same 
rate  as  Ions  are  accelerated  rearward. 

3.  Electrothermal .  The  electrothermal  thruster  Is  the 
thruster  which  fueled  this  study.  In  this  thruster,  electric 
power  or  solar  radiation  is  used  to  heat  the  propellant  to  a 
high  temperature.  The  heating  may  be  done  by  producing  an 
electric  discharge  through  the  propellant  gas  (arcjet)  or  by 
flowing  the  propellant  gas  over  surfaces  heated  using 
electricity  (reslstojet)  or  reflected  sunlight.  If  reflected 
sunlight  is  used  to  heat  the  propellant,  then  this  would 
result  in  discontinuous  low  thrust.  Discontinuous  because 
the  spacecraft  might  enter  the  earth's  shadow  during  an  orbit 
(3:1-9,  7:xiii-xvi). 

A  discontinuous,  low-thrust  orbit  transfer  between  two 
non-coplanar ,  circular  orbits  may  soon  be  accomplished  by  a 
spacecraft  currently  being  developed  at  the  Air  Force 
Astronautics  Laboratory.  The  spacecraft  is  a  solar-powered 
rocket  capable  of  transferring  large  payloads  from  low  earth 
orbit  to  geosynchronous  orbit.  The  booster  will  likely 
employ  two  solar  collectors  which  concentrate  reflected 
sunlight  and  focus  it  to  an  area  on  the  engine.  The  hydrogen 
propellant  is  heated  and  then  expelled  through  a  nozzle. 
Using  solar-heated  liquid  hydrogen  as  the  propellant,  the  low 
thrust,  high  specific  Impulse  engine  may  evolve  into  an 
efficient  powerplant  for  space  transportation  by  the  end  of 
the  century  (11:119-120). 

The  concept  of  a  solar  powered  boost  vehicle  is  not  new. 
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The  current  Air  Force  development  program  Is  an  outgrowth  of 
solar  engine  research  that  originated  in  the  early  1960's 
(11:119).  The  revival  of  this  concept  is  brought  about  by 
new  materials  and  new  control  techniques  for  the  solar 
collectors.  A  solar  powered  booster  could  be  an  important 
component  in  the  USAF  space  inventory. 

Previous  work  on  low  thrust  transfers  between  non- 
coplanar  circular  orbits  revealed  an  outward  spiral  is  the 
optimal  trajectory  (1,6:1,  10:150-154).  Cass  later  derived  a 
control  law  for  a  discontinuous,  low  thrust  transfer  between 
non-coplanar  circular  orbits  (6:16).  That  is,  a  constraint 
was  added  such  that  the  eccentricity  equals  zero.  This  study 
was  conducted  to  locate  the  optimal  launch  time  by  minimizing 
transfer  time  using  the  control  law  derived  by  Cass. 

The  derivation  is  divided  into  two  problems  of  differing 
timescales.  The  fast  timescale  problem  optimizes  the  changes 
in  the  semi-major  axis  and  inclination  over  one  orbit  with 
vehicle  mass  held  constant.  The  result  is  an  optimal  control 
law  for  one  orbit.  The  slow  timescale  problem  uses  the 
optimal  control  law  to  determine  how  much  to  change  the 
semi-major  axis  and  Inclination  on  each  orbit  so  that  the 
final  boundary  conditions  are  met  in  minimum  time.  The  mass 
and  acceleration  are  updated  after  each  orbit. 
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LLl.  derivation 


A  control  law  was  previously  derived  by  Cass  for  the 
discontinuous  low  thrust  orbit  transfer  problem(6).  The 
purpose  of  this  study  is  to  use  the  control  law  to  find  the 
optimal  launch  time.  The  optimal  launch  time  minimizes 
transfer  time.  To  avoid  duplication  of  effort,  Cass's  work 
will  not  be  rederived.  However,  the  pertinent  concepts, 
equations,  and  optimization  algorithm  will  be  fully 
described . 

The  derivation  was  divided  into  two  problems  of  differing 
timescales;  the  fast  timescale  and  the  slow  timescale.  The 
fast  timescale  problem  optimizes  the  changes  in  semi-major 
axis  and  inclination  during  one  orbit  of  the  central  body 
(earth)  with  vehicle  mass  held  constant.  The  slow  timescale 
problem  incorporates  the  fast  timescales  results  and 


determines  how 

much 

to  change 

the 

semi -major 

axis  and 

inclination  on 

each 

orbit  so 

that 

the  final 

boundary 

conditions  are  reached  in  minimum  time.  Vehicle  mass  was 
recomputed  after  each  orbit  to  account  for  propellant  loss. 

FAST  TIMESCALE  PROBLEM  (6:3-18) 

The  solution  of  the  fast  timescale  problem  will  yield  the 
optimum  thrust  control  law  for  one  orbit.  The  thrust  control 
law  maximizes  a  change  in  either  inclination  or  semi-major 
axis  when  a  particular  change  in  the  other  is  specified.  The 
derivation  assumes  two-body  motion  with  the  earth  as  the 
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central  body.  Also,  since  electric  engines  produce  low 
thrust,  the  change  in  orbital  elements  for  one  orbit  is 
assumed  to  be  very  small  (i.e.,  da  %  Aa )  .  Therefore,  the 
orbital  elements  are  considered  constant  during  each  orbit. 
The  spacecraft  mass  is  also  considered  constant  for  one  orbit 
since  propellant  loss  is  very  small. 

Lagrange's  planetary  equations  in  their  acceleration 
component  form  are  used  to  find  the  optimal  thrust  control 
law.  The  planetary  equations  are  (12:35): 

dO  W  ( 1-e2  ) 1/2  s i n  ( to* f  ) 

dt  n  a  sln(i)  (1  +  e  cos(f))  ^ 

di  W  ( 1-e2  ) 1/2  cos  ( <*>+  f  ) 

dt  '  n  a  ( 1  +  e  cos ( f ) )  ( 2  * 


dw  -u  ( 1-e2 ) 1/2  cos(f)  _  W  ( 1-e* ) 1/2  cot(i)  sin(<*>  +  f) 
dt"  nae  na(l+e  cos(f)) 

+  V  (1-e2 ) 1/2  sin(f)  (2  +  e  cos(f)) 

nae  ( 1  +  e  cos ( f  ) )  ^ 

de  U  ( 1-e2 ) 1X2  sin(f) 

dt  n  a 


+ 


V  (1-e2 ) 1/2 
nae 


1  +  e  cos(f) 


(1-e2) 
1+e  cos(f) 


(4) 
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da 

dt 


2  U  e  sin( £  ) 


+ 


2  V  ( 1  +  e  cos ( £  )  ) 


(5) 


n 


. ,  2.1/2 
( 1-e  ) 


n 


. ,  2 . 1/2 
( 1-e  ) 


where  £2  is  the  longitude  of  the  ascending  node,  w  is  the 
argument  of  perigee,  i  is  the  spacecraft  inclination,  e  is 
the  orbit  eccentricity,  a  is  the  semi-major  axis,  and  £  is 
the  true  anomaly.  These  orbital  elements  are  shown  in  Figure 
1.  The  variables  U,V,  and  W  are  the  radial,  tangential,  and 
normal  acceleration  component,  respectively,  and  n  is  the 
mean  motion  which  is  given  by: 

n  *  ( M/a*) 1/2  (6) 

where 

M  =  G  (7) 

and  G  is  the  universal  gravitation  constant  and  is  the 
mass  of  the  earth. 

For  this  problem,  the  initial  and  final  orbits  are 
defined  to  be  circular  (e  =  0).  In  addition  to  these 
boundary  conditions,  the  eccentricity  for  each  orbit  is  set 
equal  to  zero  throughout  the  entire  orbit  transfer.  That  is, 
a  constraint  is  introduced  which  requires  the  change  in 
eccentricity  between  each  orbit  be  equal  to  zero.  This 
constraint  greatly  simplifies  the  planetary  equations,  as 
shown  below: 


dft  W  sin(w+f) 

dt  ~  n  a  sln(i) 
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di 

dt 


W  cos ( W+f  ) 
n  a 


(9) 


dw 

—  =  undefined  for  a  circular  orbit  (10) 


de  U  sin(f)  +  2  V  cos(f) 

dt  n  a  n  a 


(11) 


da  2  V 

dt  n 


(12) 


In  general,  w  is  undefined  for  a  circular  orbit.  Thus, 
the  new  reference  point  for  f  is  the  ascending  node  as  shown 
in  Figure  1.  For  continuous  thrust  vehicles,  f  is  the 
vehicle's  angular  displacement  with  the  engines  on.  In  this 
problem,  thrust  is  posjible  only  in  the  sunlight. 
Accordingly,  angle  f  is  redefined  as  the  angular  displacement 
of  the  spacecraft  in  the  sunlight  measured  with  respect  to 
the  earth's  shadow  exit  point.  A  new  angle,  s,  is  introduced 
to  give  the  angular  displacement  between  the  ascending  node 
and  the  shadow  exit  point.  Now,  there  is  a  complete 
definition  of  the  spacecraft's  position  with  respect  to  the 
ascending  node  as  shown  in  Figure  2.  Thus,  the  angle  (w+f) 
in  Eqs  (8)  and  (9)  is  replaced  by  (s  +  f).  Also,  the  angle  g 
is  measured  between  the  shadow  exit  to  the  next  shadow  entry 
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Figure  2.  Shadow-Related  Angles 


10 


and  m  is  one-half  of  the  shadow  angle.  See  Figure  2  for  a 
depiction  of  these  angles.  Uses  of  angles  g  and  m  are 
discussed  later  in  this  section. 

Additional  modifications  to  the  planetary  equations 
include  changing  the  independent  variable  from  t  to  f  (using 
n  =  df/dt)  and  replacing  n  with  (p/a9)1/2  (Eq  (6)).  The 
planetary  equations  are  now: 


dn  W  a2  sin(f+s) 

df  ”  p  s in ( 1 ) 


(13) 


di  W  a2  cos(f+s) 

df  =  p 


(14) 


de  a2 

__  _  — —  [u  s in ( f )  +  2  V  cos ( f )  I  (15) 

dr  p 


da  2  V  a9 

df  =  p 


(16) 


To  determine  the  changes  in  the  orbital  elements  for  one 
orbit,  Eqs  (13)  through  (16)  should  be  integrated  from  0  to 
2n;  but  U,  V,  and  W  are  zero  when  the  vehicle  is  in  the 
earth's  shadow.  Therefore,  these  equations  are  integrated 
from  0  to  g,  where  g  is  measured  from  the  shadow  exit  to  the 
next  shadow  entry. 


g  *  2n  -  2m  *  2(rr  -  m) 


(17) 


where  m  is  one-half  of  the  shadow  angle(see  Figure  2). 
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The  accelerations  u,  V,  and  w  are  modeled  as  functions  of 


f .  From  Figure  3 : 


U  =  T  cos y  cos a 

V  =  T  cosy  sina 

V  =  T  siny 


(18) 
(19  ) 
(20) 


where 

y  »  y(£)  (21) 

a  =  a(f)  (22) 

and  T  Is  the  total  acceleration  due  to  the  thrust. 
Substituting  Eqs  (18)  through  (20)  Into  Eqs  (13)  through  (16) 
and  integrating  between  0  and  g,  the  change  In  orbital 
elements  become: 

9 

(23) 


T  a  r 

AO  *  dn  =  ■ ^  sin( t )  I  sin(y)  sin(f+s)  df 


T  a 


Ai  %  di  = 


Jsin 


(y)  cos(f+s)  df 


(24) 


T  a 


Ae  %  de  s 


i 


(cos (y)  cos(a)  sin(f) 


+  2  cos(y)  sln(a)  cos(f))  df 


(25) 
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Aa  %  da 


cos (r)  sin(a)  df 


(26) 


2  T  a 


J 


The  following  conditions  are  used  to  find  the  optimal 
thrust  control  law  for  one  orbit:  di  is  maximized  for  a 
given  da,  subject  to  the  constraint  de  =  0 .  The  performance 
index  with  constraint  relationships  is: 


J  ( a,  y )  =  j 


T  a 


sin(r)  cos(f+s)  df 


+  X 


I 


2  T  a 


cos (r)  sin(a)  df 


da 


L  o 


+  X 


J 


T  a‘ 


*  J  V 

O 


(cos (r)  cos(a)  sin(f)  +  2cos(y)  sin(o)  cos(f))  df 

(27) 


where  X^  and  X^  are  Lagrange  multipliers.  The  goal  is  to 
find  the  functions  y  and  a  that  produce  a  stationary  value  of 
J.  The  solution  of  this  problem  was  done  by  Cass  and  is 
patterned  after  an  optimization  scheme  given  by  Bryson  and  Ho 
(4:48-49).  The  optimal  thrust  control  law  is: 


a  =  tan 


2 (X^a  +  X^cos ( f ) ) 


X^  sin(f) 


(28) 


tan 


-i 


cos (f +s ) 


fctX*  sin* ( f )  +  4(Xta  +  Xjcos(f ) )*]*'* 


(29) 
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Substituting  Eqs  (28)  and  (29)  Into  Eqs  ( 24)  through  (26) 
yields : 


T  a‘ 


Ai  = 


l 


cos  (f+s)  df 


p  J  (X*  sin*(f)  +  4(X  a  +  X  cos(f))*  +  cos2(f+s)J 
0  2  12 


1/2 


(30) 


_  z  9  (X  sin  (f)  +  4cos(f)  (X  a  +  X  cos(f))  df 
i  a  r  z  12 


Ae 


J 


M  J  [X2  sin2 ( £ )  +  4 (X  a  +  X  cos(f))2  +  cos2 ( f+s )] 1/2 
0  2  12 

(31) 


4  T  a 


Aa  = 


I 


[X^a  +  X^  cos ( f )  ]  df 


^  J  [X2  sin2 ( £ )  +  4 (X  a  +  X  cos(f))2  +  cos2 ( f+s )) 1X2 
0  2  12 

(32) 


Assuming  a  sperical  earth.  Ail  is  neglected  and  will  not  be 
considered  in  this  problem.  For  convenience,  let: 

u  =  X^a  (33) 

and 


D  =  X2  sin2(f)  +  4(u  +  Xcos(f))2  +  cos2(f+s).  (34) 
2  2 


Substitution  yields: 
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*  „  s  9  lu  +  A  cos(f )  ]  df 

4  T  a  r  2 

a  =  -  - —  (37) 

"  J  D1/2 

O 

Eqs  (35)  through  (37)  represent  changes  in  i,  e,  and  a 
for  one  orbit  using  the  optimal  thrust  control  law. 
Actually,  this  control  law  will  maximize  Ai  for  a  given  Aa, 
subject  to  the  constraint  Ae  =  0 .  Eqs  (35)  and  (37)  will  be 
used  in  the  slow  timescale  problem  to  ensure  optimal  control 
during  each  orbit  and  to  define  the  amount  of  change  possible 
for  Ai  and  Aa  on  a  given  orbit.  The  Ae=0  constraint  is 
solved  numerically  to  determine  u  and  \^( further  use  of  Eq 
(36)  is  given  in  the  next  section). 

However,  before  these  equations  are  implemented  the 
angles  m  and  s  must  be  calculated.  Recall  m  is  one-half  of 
the  shadow  angle  and  s  is  an  angle  in  the  orbit  plane  between 
the  ascending  node  and  the  point  where  the  spacecraft  exits 
the  earth's  shadow.  Sun-earth  geometrical  relationships  are 
used  to  determine  m  and  s.  According  to  Link  (9:122-127), 
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the  following  sun-earth  geometrical  relationships  are  used  to 
calculate  m  (see  Figure  4): 


cos ( m) 


cos(a  ) 

O 

COS (B ) 


(38) 


where  B  and  a  are  angles  shown  in  Figures  4  and  5  and  a  is 

©  O 

half  of  the  maximum  shadow  possible  for  a  given  altitude. 

Link  also  gives  expressions  for  the  angles  B  and  a  : 

© 

T  4/? 

cos(B)  =  II  -  (cos  (<50)sin(  i  )sin(aQ-n)  -  sin  (<50)cos  (  i  ) )  J 

(39) 

and 

cr  =  n  +  n  -  R  (40) 

o  O  ■  O 

or 

cos  (a  )  =  F  II  -  a"*!4'*  ♦  IG  /a)  (41) 

o  o  o 

where  o»Qand  <50  are  the  right  ascension  and  declination  of  the 
sun  in  earth  centered  coordinates  and  F  and  G  are  constants 

o  o 

depending  on  the  size  of  the  earth  and  sun  and  the  distance 
between  them. 

F  =  0.99994668  G  =  0.0046229536 

O  O 

Link  also  has  an  expression  for  s: 


where 


s  =  (tt/2)  -  L  +  m 

cos (6g) 

■  -cos(a0-0) 


(42) 


L 


cos(B) 


(43) 


All  ingredients  should  now  be  available  for  the  slow 
timescale  problem. 


SLOW  TIMESCALE  PROBLEM  (6:26-37) 


The  purpose  of  solving  the  slow  timescale  problem  is  to 
determine  how  much  to  change  the  semi-major  axis  and 
inclination  on  each  orbit  such  that  final  boundary  conditions 
are  met  in  minimum  time.  Fast  timescale  results  for  Aa,  Ai, 
and  Ae  are  used  to  ensure  optimal  control  during  each  orbit 
and  to  define  the  amount  of  change  possible  for  Aa  and  Ai  on 
a  given  orbit.  Mass  is  recomputed  for  each  orbit  to 
compensate  for  propellant  loss. 

Before  the  minimum  time  control  problem  is  solved,  a  few 
preliminary  steps  are  taken.  An  expression  is  needed  for 
da/dt  and  di/dt  for  the  slow  timescale  problem.  Since  the 
orbital  elements  change  very  little  on  each  orbit,  these 
rates  are  approximated  by: 


da  Aa 

dt  *  At 


(44) 


di  Ai 

dt  *  At 


(45) 


where  At  is  the  orbital  period.  For  a  circular  orbit. 

At  «=  2tt  la'/nl*'*-  (46) 


However,  the  elements  a  and  i  only  change  when  the  spacecraft 
is  in  sunlight.  Therefore,  At  must  be  altered  to  only  cover 


20 


that  portion  of  the  orbit.  The  new  At  is: 

^  S  l/i 

At  *  ~ 2n~  period)  =  2(n-m)[a  /v 1 

Dividing  Eqs  (35)  and  (37)  by  At  yields: 


da 

dt 


2  T  a 


i/i  .  . 

V  (n-m) 


(u  ♦  X^cosff)  df 

V7* 


(47) 


(48) 


di  T  a1/l  r  cos*(f+s)  df 

dt  2m1/2  (n-m)  oJ  D4'* 


(49) 


Since  the  goal  is  to  find  a  minimum  thrusting  time,  only  the 
time  in  which  the  spacecraft  is  in  the  sunlight  is 
considered.  The  total  thrusting  time  for  the  entire  orbit 
transfer  is  the  summation  of  thrusting  time  for  each  orbit. 

As  previously  stated,  the  spacecraft's  mass  is  recomputed 
for  each  orbit  to  account  for  propellant  loss.  The 
propellant  loss  causes  an  increase  in  acceleration.  For  a 
constant  thrust  ion  rocket,  propellant  flow  is  constant. 
Acceleration  can  be  modeled  as  a  function  of  time: 


I 

e 

T ( t )  *  - : -  (50) 

1  -  m  t 

p 


where  T  is  the  initial  spacecraft  acceleration,  t  is  the 

O 

time,  and  m  is  the  specific  mass  flow  rate  (mass  flow  rate 

P 
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divided  by  Initial  spacecraft  mass). 

Eqs  (48)  and  (49)  can  be  solved  more  easily  If  there  is 
no  time  dependence.  To  eliminate  the  time  dependence,  a 
change  of  variable  is  made.  The  new  independent  variable  t 
is  defined  by; 


where  t  is  the  total  accumulated  velocity  change.  The  units 
of  r  are  DU/TU.  Minimizing  t  will  minimize  the  thrusting 
time  and  fuel  expended.  Using  t  as  the  independent  variable 
will  not  yield  the  number  of  orbits  and  total  thrusting  time 
(in  seconds)  when  the  final  boundary  conditions  are  met. 
These  values  depend  on  the  engine  performance  parameters  (T  , 

O 

• 

mp  ,etc...).  The  solution  using  t  as  the  independent 
variable  will  give  the  change  in  velocity  needed  to  reach  the 
final  boundary  conditions  for  any  type  of  low  thrust  engine. 
By  integrating  Eq  (51),  an  expression  relating  the  engine 
performance  parameters  and  time  to  t  is  found.  With  r  =  0 
when  t  =  0 

•  • 

t  *  -(T  /  m  )  In  (1  -  m  t)  (52) 

op  p 

Substituting  Eq  (51)  into  Eqs  (48)  and  (49)  gives: 
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da 

dr 


di 

dr 


2 

1/2 

V 


S/2 

a 


(n-m) 


1/2 

a 


2fj*'x  (n-m) 


9 


f 


9 


1 


o 


(u  +  X^cos ( f ) )  df 


D 


1/2 


cos2(f+s)  df 

V7* 


(53) 


(54) 


Eqs  (53)  and  (54)  are  essential  in  finding  the  minimum 
time  solution.  Cass  provides  the  derivation  to  the  minimum 
time  problem  using  the  minimum  time  solution  given  in  Bryson 
and  Ho  ( 4 : 87-89 )  . 

Bryson  and  Ho  have  shown  a  minimum  time  solution 
satisfies  the  following  three  conditions: 

H(tf)  »  0  (55) 

K  =  1/2,3, . . .  ,r  (56) 

j  =  1,2,3, ...,q  (57) 

where  H  is  the  Hamiltonian,  X^'s  are  Lagrange  multipliers, 
x.'s  are  the  state  variables,  and  the  u.  's  are  the  control 
variables.  The  prime  indicates  differentiation  with  respect 
to  t.  In  this  problem,  r  *  1  (u^  *  u)  and  q  =  2  (x^  *  a  and 
x^  *  1).  Eqs  (55)  through  (57)  will  now  be  applied  to  this 
specific  problem. 


d  H 

9  u. 


=  0 


X  = 


-  9  H 

9  x 
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Foe  this  problem,  the  Hamiltonian  is 


da  di 

H  =  \  — —  +  V  —2 —  +  1 
*  dT  v  dT 


(58) 


Since  t  does  not  appear  explicitly  in  H, 


H  =  0 


Therefore,  H(tf)  =  0  Implies 


H  ( t )  =  0 


(59) 


(60) 


for  all  t  i  0.  So  Eq  (58)  becomes 


a  dT 


di 

+  \  -ST-  +  1 


(61) 


for  all  time.  Substituting  Eqs  (53)  and  (54)  into  (61) 


gives : 


j-  2  aS/S 
[  S'*  ( n-i 


-[ 


(u  +  \  cos(f)  df 

Z 


lS'*  (n-m) 


9  _ 

r  cos*(f+s)  df  I 

j  — ^ —  r  l-° 


(62) 


For  this  problem,  Eq  (56)  becomes 


♦  v  f  ii  I 

u  dr  0  u  L  dT  J 


(63) 
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Cass  derives  the  contents  of  Eqs  (68)  and  (69).  Substituting 
Eqs  (67)  through  (69)  into  (65)  yields 


r  a1'2  [ 3 ( n-m)  +  2  X  a) 


X  «  -X 
a  a 


(U  +  X  COS ( f ) )  d£ 


i/*  .  .2 

/u  (n-m) 


r  iu  r  » 

J  D 


1/2 


2  X  a 


»/2 


1/2  .  . 

V  (rr-m) 


_  r  *  "2 

-m)2  J 


(u  ♦  X  cos(f))  sin(f-»-s)  cos(f+s)  df 


»/2 


8/2 


4  X  a  (u  +  X  cos ( 2m) ) 


m)  (X2sin2(2m)  +  4 ( u+X#cos  (  2m) ) 2 ♦  cos*  (s-2m)  ) 1/2 


-X 


(n-m)  +  2  X  a 


4  v'*  a1'2  (n-m) 


^  J 


cos  (£+s)  df 


.1/2 


X  a 


1/2 


1/2  ,  4 

V  (n-m) 


J 


sin(f+s)  cos(f+s)  df 


.1/* 


X  a 


1/2 


2  (n-m) 


I 


sin(f+s)  cos  (f+s)  df 


.8/2 


X  a1/2  cos2(s-2m) 


Ai*'2(n-m)  lX*sin*(2m)  +  4 ( u+X^cos ( 2m) ) *+  cos*(s-2m)) 


1/2 


where 


G  -  F  (a2-l)"1/2 

o  © 


2  ,  2 .  2 .  . , 1/2 
a  (cos  (B)  -  cos  (cr  ) ) 

O 


(70) 


(71) 
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Now,  Eqs  (62)  and  (64)  can  be  solved  simultaneously  for  \ 

a 

and  V  in  terms  of  integrals  containing  the  state  and  control 
variables.  From  Eg  (64)/ 


X 


X  a 

a 


J  - 


(X*  sin* ( f )  +  cos* ( f +s ) ]  df 


»  2 
r  (U  +  X^  COS(f))  COS  (f+S)  df 

j  ^s/* 


(74) 


Substituting  this  result  into  Eg  (62)  gives 


-1 

-  x/z  .  . 

2  (j  (rr-m) 

»  2 
•  (u  +  Xcos(f))  cos  (f+s)  df 

Z 

Q(U/S) 

z/z 

a  oJ 

D  »'* 

(75) 


where 


Q( 


U/S)=J 


g  9  z 

4(u  +  X  cos(f))  df  ( u  +  X  cos ( f ) )  cos  (f+s)  df 


,t/i 


r  iu  t  ^co5  mi  ci 

J  n*"* 


9  9 

cos*(f+s)  df  r  (*.*  sin*(f)  +  cos*(f+s)l  df 


r  cos  iitsi  ar  p 

J  ~[Tz  J 


(76) 


Egs  (74)  and  (75)  can  be  substituted  into  either  Eg  (70)  or 
(72)  to  also  produce  an  integral  equation  in  terms  of  the 
state  and  control  variables.  Using  the  new  expression  for 
Egs  (70)  or  (72)/  the  two  state  Egs  (53)  and  (54),  and 
equating  Eq  (36)  to  zero  to  find  X^,  the  minimum  time  problem 
can  be  solved. 
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III.  SOLUTION  ALGORITHM 


The  solution  to  the  minimum  time  problem  uses 
dimensionless  parameters.  For  instance,  a  is  measured  in 
DU's,  r  is  measured  in  DU's  and  TU's,  the  gravitational 
parameter  fj  is  equal  to  1  DU*/TUl,  and  all  angles  are 
measured  in  radians. 

Choose  either  Eq  (70)  or  (72)  to  use  in  the  solution.  In 
this  solution,  equation  (70)  is  used.  It  might  appear  Eq 
(72)  is  less  complicated  since  it  has  fewer  terms  and  Y  is  a 
factor  in  all  terms.  Notice,  X  is  a  constant  when  m  -  0  and 

X. 

Y  =  0.  Cass  revealed  that  in  the  range  of  interest  for  X^,  u 
is  double-valued  (6:39).  It  was  also  found  X  varied 
dramatically  for  different  values  of  a  and  i.  See  Figures 
(6)  and  (7).  Notice  the  X  curve  maintains  its  shape  and 

CL 

only  changes  slightly  between  initial  and  final  conditions. 
The  X.  curve  is  not  as  well  behaved.  The  fact  u  is 

V 

double-valued  for  X^  causes  the  search  scheme  to  diverge. 
That  is,  for  a  particular  value  of  X^  ,  there  could  be  two 
values  of  u.  Which  value  of  u  is  better  and  how  can  one  use 
a  numerical  search  scheme  to  find  that  u?  To  avoid  this 
complication,  equation  (70)  is  used.  Now,  for  a  particular 
X^,  there  is  only  one  value  of  u. 
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The  minimum  time  solution  is  calculated  in  the  following 


way: 

1.  Input  initial  boundary  conditions. 

initial  a  =  1.03  DU 

initial  i  =  0.4974  rad  (28.45  degrees) 

2.  Pick  an  initial  value  for  X.  (X  (0)).  See  Figures  8 

a  a 

through  12  for  X^(0)  during  the  analyzed  dates. 

3.  Determine  sun-earth  parameters 

input:  a,  i,  fl,  day 

output :  m,  B,  L,  a .  Then  calculate  s  and  g. 

4.  Iterate  equation  (75)  in  order  to  find  what  u  would 

produce  that  X^.  See  step  4a  concerning  X^  in  equation 

(75).  NOTE:  Simpson's  rule  was  used  to  evaluate  the 

Integrals  (5:165).  Bisection  method  was  used  for  the 

iteration  of  u  (5:29). 

lUfijlt:  X  ,  s,  g,  m,  a 
<* 

autfillt:  u,  X. 

4a.  \z  was  calculated  for  each  guess  of  u  using  the 

constraint  Ae  =  0  (Eq  (36)).  The  bisection  method  was 
used  to  determine  X^. 
input :  u,  a,  g 

QUtPUt :  \ 
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Initial  Lambda 


Ascending  Node  (degrees) 


\q(0)  During  the  First  Day  of  Winter 


X^(0)  Midway  Between  Winter  and  Spring 


Figure 


Figure 


10.  X^(0)  During  the  First  Day  of  Spring 


11. 


X 


(0)  Midway  Between  Spring  and  Summer 
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5.  Since  small  changes  in  orbital  elements: 
new  a  *  a  +  Aa 
new  1  *  i  -  Ai 

where  Aa  and  Ai  are  calculated  using  Eqs  (53)  and  (54), 
respectively.  These  equations  were  not  integrated  since 
the  changes  are  very  small.  To  approximate  the 
integration.  At  was  £ixed  at  0.01. 
input :  u,  \2,  a,  m,  s,  g 
output :  A a,  Ai 

6.  Use  Eq  (70)  to  £ind  the  new  X 

a 

new  \  *  X  ♦  AX 

a  a  a 

where 

,  dX  AX 

X^  a  — a  9i  *  rhs(right  hand  side)  of  Eq  (70) 

then 

AX  =  (rhs)  At 

a 

The  rhs  of  Eq  (70)  uses  the  current  value  of  a  (not  the 
new  a  calculated  in  step  5).  At  «  0.01. 

in PUl:  X  ,  X  ,  a,  m,  s,  g,  X  ,  G  ,  P  ,  B,  a  ,  u 

O  l  z  o  o  o 

output:  new  X^ 

7.  Check  if  met  final  boundary  conditions  where: 

final  a  «  6.6  DU 
final  i  ■  0.0  radians 
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Use  the  flowchart  for  checking  final  boundary  conditions 


where  tol  ■  0.001  radian 
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1Y.  RESULTS  AHJ& 


This  study  analyzes  the  transfer  times  for  key  dates 
during  one-half  of  a  year  from  the  first  day  of  winter  to  the 
first  day  of  summer.  The  transfer  times  for  the  other  half 
of  the  year  are  assumed  to  be  similar  to  those  of  the 
analyzed  dates.  The  analyzed  dates  are:  first  day  of 

winter (21  DEC)/  midway  between  winter  and  sprlng(4  FEB)/ 
first  day  of  spring (21  MAR)/  midway  between  spring  and 
summer (6  MAY)/  and  the  first  day  of  summer (21  JUN).  The 
reference  day  used  in  the  computer  program  is  the  vernal 

equinox,  resulting  in  the  analyzed  day  numbers  to  be  274, 

320/  0/  46/  and  91,  respectively.  A  point  should  be  made 

concerning  the  number  of  data  points.  The  number  of  data 

points  per  day  is  chosen  to  yield  general  characteristics  of 

the  data  and  is  not  designed  to  give  exact  numerical 

quantities . 

All  transfers  required  the  same  change  in  the  semi-major 
axis  and  inclination.  The  Initial  bounday  conditions  are  a ■ 
200  km  *  1.03  DU  and  i.=  26.5  degrees  =  0.4974  radians  (a 
possible  shuttle  orbit).  The  final  boundary  conditions  are 
af  =  6.6  DU  and  if  =  0.0  radians.  Recall  the  total 

accumulated  velocity  change  (t)  is  minimized  rather  than  the 
actual  transfer  time.  If  one  minimizes  t  then  transfer  time 
is  also  minimized. 

Figures  13  through  17  show  how  t  varied  with  n  during 
each  of  the  five  analyzed  days.  Each  plot  of  t  Is  a 
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0  BO  180  270  300 

Ascending  Node  (degrees) 


Figure  13. 


t  on  the  First  Day  of  Winter 


Figure  14.  t  Midway  Between  Winter  and  Spring 
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0  90  180  270  300 

Ascending  Node  (degrees) 

Figure  15.  t  on  the  First  Day  of  Spring 


Figure  16.  t  Midway  Between  Spring  and  Summer 


"sinusoid*'  of  sorts  and  varies  in  a  predictable  manner.  The 
plot  for  spring  (Figure  15)  is  symmetric  for  t.  The  minimum 
r  on  each  day  in  Figures  13  through  17  is  plotted  versus  the 
day  in  Figure  18.  The  ascending  node  at  which  the  minimum 
t's  occurred  are  also  plotted  versus  the  day  in  Figure  19. 
The  following  paragraphs  will  examine  Figures  13  through  19. 

An  important  trend  arises  concerning  the  optimal 
ascending  node  in  Figures  13  through  17.  Note:  the  optimal 
ascending  node  is  defined  as  the  ascending  node  where  t  is  a 
minimum.  The  line  of  nodes  (a  line  connecting  the  ascending 
and  descending  node)  for  the  minimum  t  orbit  on  each  day  is 
perpendicular  to,  or  nearly  perpendicular  to,  the  earth-sun 
line.  A  comparison  of  Figures  13  through  17  and  20  through 
24  Indicates  the  optimal  launch  time  (i.e.,  acsendlng  node) 
is  the  one  which  minimizes  the  time  in  the  earth's  shadow. 
Figures  20  through  24  show  the  ascending  node  for  the 
quickest  escape  of  the  earth's  shadow  occurred  at,  or  very 
near,  the  ascending  node  for  the  minimum  t  orbit.  The 
ascending  node  for  the  minimum  initial  shadow  half-angle  (m  ) 

O 

also  corresponds  to  the  minimum  t  orbit.  Obviously,  the  less 
time  spent  in  the  earth's  shadow  allows  more  thrusting  time 
per  orbit.  More  thrusting  time  per  orbit  yields  smaller 
transfer  times. 

Assuming  the  minimum  t  orbit  requires  the  line  of  node  to 
be  perpendicular  to  the  earth-sun  line,  the  question  remains 
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Figure  21.  Effects  of  (1  on  m  Midway  Between 
Winter  and  Spring 
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as  to  which  node  Is  the  ascending  node.  If  the  northern 
hemisphere  is  between  fall  and  spring,  then  the  ascending 


node  is  (aQ  -  n/2).  If  the  northern  hemisphere  is  between 
spring  and  fall,  then  the  ascending  node  is  (aQ  +  n/2).  For 
the  first  day  of  spring  and  fall,  the  optimal  ascending  node 
is  both  (ct  ±  n/2).  These  orientations  will  minimize  m  , 

Q  o 

thereby  minimizing  t.  Figure  19  shows  how  the  ascending  node 
at  minimum  t  changed  with  respect  to  the  day  of  the  year. 
The  right  ascension  of  the  sun  is  also  included  in  Figure  19 
for  comparison  purposes.  The  two  optimal  ascending  nodes  for 
spring  is  due  to  the  symmetric  lighting  conditions. 

Another  trend  in  the  data  is  the  change  in  semi-major 
axis  (Aa)  is  greater  than  the  change  in  inclination  (Ai) 
during  the  initial  phase  of  the  transfer.  Figure  25  displays 
a  typical  profile  for  Aa  and  Ai.  This  agrees  with  the  fact 
that  plane  changes  require  much  more  energy  than  coplanar 
changes  when  gravity  is  a  major  factor  (1,  2:169-170).  This 
profile  is  also  the  path  needed  by  the  spacecraft  to  minimize 
or  eliminate  (depending  on  the  time  of  year)  the  time  in  the 
earth's  shadow. 

Figure  18  shows  the  overall  minimum  r  orbit  occurring  at 

both  the  first  day  of  summer  and  the  first  day  of  winter. 

This  corresponds  with  the  overall  minimum  m^  which  is  plotted 

in  Figure  26.  The  difference  between  the  t's  in  Figure  6  is 

only  1.21%.  Whereas  the  change  in  n^'s  in  Figure  26  is 

10.91%.  However,  these  differences  in  t  and  m  are  not 

© 
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Figure  25. 


Propagation  of  a  and  i  During  an  Orbit 
on  4  FEB  88  with  Q  =  180  degrees 
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a  o  (DU) 


consistent  on  any  particular  day  (Figures  13  through  17  and 
20  through  24).  For  instance,  t  varied  from  4.7%  to  5.9%  on 


any  of  the  analyzed  dates.  Whereas  »o  varied  from  2.6%  to 
13.8%.  The  upper  limit  of  these  changes  for  r  and  m^ 
occurred  simultaneously;  as  did  the  lover  limit.  Table  I 
contains  the  minimum  and  maximum  values  of  t  and  m  for  the 

O 

analyzed  dates. 


TABLE  I.  Changes  in  t  and  ra^ 


Date 

(day#) 

min  1 
max  J 

(DU/TU) 

min  \ 
max  J 

(radians ) 

%  change 
in  t 

%  change 
in  m 

O 

21  DEC 
(274) 

0.77232 

0.81745 

1.163663 

1.324214 

5.80 

13.80 

4  FEB 

(320) 

0.77865 

0.81628 

1.222310 

1.324079 

4.83 

8.32 

21  MAR 
(0) 

0.78165 

0.81850 

1.290653 

1.324239 

4.71 

2.60 

6  MAY 
(42) 

1.220574 

1.324144 

5.07 

8.49 

21  JUN 
(91) 

wm 

1.163673 

1.324218 

5.90 

13.80 
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CQMCLU&IQHS  AttE.  RECOMMENDATIONS 


This  study  determined  the  optimal  launch  time  for  a 
discontinuous,  low  thrust  orbit  transfer.  The  optimal  launch 
times  are  the  first  day  of  winter  and  summer  and  occurred 
simultaneously  with  minimum  m  (smallest  initial  shadow  half- 

O 

angle).  The  actual  transfer  time  in  seconds  is  not 
calculated  because  the  independent  variable  was  changed  from 
t  to  r  (the  accumulated  velocity  change).  This  was  done  in 
order  to  analyze  the  entire  class  of  low  thrust  rockets.  If 
one  wants  to  calculate  the  transfer  time  in  seconds  for  a 
specific  rocket,  then  use  Eq  (51)  or  (52)  to  change  the 
independent  variable  from  t  to  t.  If  the  independent 
variable  is  t,  then  the  changes  in  aQ  and  6Q  during  the 
transfer  should  be  taken  into  account. 

The  following  is  a  rough  calculation  for  converting  from 
r  to  t.  Solve  for  t  using  Eq  (52) 


where 


t 


exp(-T  m  /  T 

P  O 


« 


thrust  force  F 

To  *  vehicle  mass  *  m  . 

v«h 


m 


propellant  mass  flow 
vehicle  mass 


v  •  h 
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As  an  example,  use  the  engine  parameters  given  in  Re£  11: 

F  *  0.8  lb 

m  =  0.001  lb/sec 

with  a  payload  o£  10,000  lb  (a  large  payload).  Assume  the 
total  vehicle  weight  (including  payload)  is  14,000  lb.  Then 

F  0.8  lb 

T  =  - -  =  - -  *  1 . 84x10  *  ft/sec* 

°  v«h  14000  lb  /  32.2  ft/sec 


.  m  0.001  lb/sec 

mp  =  m  :  =  14000  lb 

v*n 


7.14  x  10‘*  /sec 


r  .  =  0.77232  DU/ TU  * 

m\.  n 


2.093X107  ft 
1  OU 


1  TU 

*  - 

806.8  sec 


=  20,035  ft/sec 


Substituting  yields  t  «  7,566,926  sec  «  2102  hr  =  88  days. 
Using  t  of  0.8185  DU/TU  yields  t  *  91  days;  a  3.9% 

WMM 

increase  in  transfer  time. 

As  a  continuation  to  this  study,  one  could  find  the 
optimal  launch  time  without  constraining  Ae  to  be  equal  to 
zero.  Another  topic  of  interest  is  to  take  into  account  a 
non-spherical  earth,  thereby  causing  a  regression  in  the  line 
of  nodes  and  rotation  of  the  line  of  apsides  (if  orbit  is 
eccentric) . 
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